---
title: "Delegative Peacebuilding"
author: "Sally Sharif"
date: "2023-12-11"
output: pdf_document
---

This R Markdown file contains all replication code for "Delegative Peacebuilding" (2025). The code is organized into labeled sections for ease of navigation; for example, the first section is titled {r ImportData} and the second {r MatchPropensityScores}. The sections contain the code used to generate the tables and figures presented in both the main manuscript and the online appendix. The also present results obtained with alternative R packages.

```{r ImportData}
library(readxl)
library(tidyr)
library(dplyr)
library(arm)
library(optmatch)
library(RItools)
library(MatchIt)
library(table1)
library(corrtable)
library(glmmTMB)
library(vcd)
library(MASS)
library(DHARMa)
library(lmtest)
library(car)
library(texreg)
library(HonestDiD)
library(fixest)
library(sandwich)
library(colmaps)
library(sf)
library(ggplot2)
library(rgdal)
library(rgeos)
library(scales)
library(gridExtra)
library(performance)

data0 <- get(load("assassinations_2014_2020.RData"))
nrow(data0) # should be 1379

data <- get(load("municipal_data.RData"))
nrow(data) # Should be 1,121
```

```{r AssassinationData}
# getting number of leaders killed monthly
data2 <- data0 %>%
  group_by(codigo, year, month) %>%
  summarize(leader_killed = n()) %>% ungroup()
# add all post-conflict months to the data and fill with 0 if no value
# 0 means no leader was killed in a given municipality in a given month
data3 <- data2 %>%
  complete(codigo, year = 2014:2020, month = 1:12, fill = list(leader_killed = 0))
# getting number of leaders killed annually
data4 <- data0 %>%
  group_by(codigo, year) %>%
  summarize(leader_killed = n()) %>% ungroup()
# add all post-conflict years to the data and fill with 0 if no value
# 0 means no leader was killed in a given municipality in a given year
data5 <- data4 %>% 
  complete(codigo, year = 2014:2020, fill = list(leader_killed = 0)) 
```

```{r MatchPropensityScores}
# matching on propensity scores
# using the log because coca (log) is used as a control variable in the main analysis
data$log_coca2016 <- log(data$coca2016 + 1)
# for reproducibility
set.seed(123)
# fit the bayesian logit model with pdet as the outcome
bayesglm <- bayesglm(pdet ~ poverty + log_coca2016 + inst_development + inst_performance, data=data, family = binomial(link = "logit"))
summary(bayesglm)
## exact match
psdist <- match_on(bayesglm, data = data)
dist2 <- psdist + exactMatch(pdet ~ conflict_intensity, data = data)
# full match
fm0 <- fullmatch(dist2, data = data) 
effectiveSampleSize(fm0)
summary(fm0, min.controls = 0, max.controls = Inf, propensity.model = bayesglm)
## matched table
# Appendix Table 16
summary(dist2)
ftable(Class = data$conflict_intensity, Trt = data$pdet, fm0, col.vars = c("Class", "Trt"))
# balance table
# Appendix Table 15 and 17
xbtotal <- xBalance(pdet ~ poverty + log_coca2016 + inst_development + inst_performance + conflict_intensity, strata = list(unstratified=NULL, full=~fm0), data = data, report = "all")
print(xbtotal)
# Appendix Figure 7
png("Balance table.png", width = 4000, height = 3000, res = 300)
plot(xbtotal, variable.labels = c(poverty = "poverty", log_coca2016 = "coca (log)", inst_development = "institutional development", inst_performance = "institutional performance", conflict_intensity = "conflict intensity"))
dev.off()
# extract matched data
data$fm0 <- fm0
```

```{r MatchCEMcriminal}
## coarsened exact matching for criminal pathway variables (Albarracin et al. 2023)
set.seed(123)  
cem_match <- matchit(
  pdet ~ log_coca2016 + conflict_intensity, 
  data = data, 
  method = "cem")
# summarize the matching results
# Appendix Tables 19 and 20
summary(cem_match)
# extract balance data 
balance_data <- summary(cem_match)$sum.matched
# extract matched data
data$cem_subclass <- cem_match$subclass
```

```{r MatchCEMpolitical}
# Coarsened exact matching with Albarracin et al. (2023) variables
# exponentiate logleft2015 to remove the log for matching
data$left2015 <- exp(data$logleft2015)
# remove missing values: some municipalities in Amazonas and Vaupes have missing values.
cleaned_data <- data[complete.cases(data$logleft2015, data$participation, data$nep_mun), ]
# check the number of rows before and after cleaning
nrow(data)  # Original number of rows
nrow(cleaned_data) # Number of rows after removing missing data
# perform Coarsened Exact Matching
set.seed(123)  
cem_match2 <- matchit(
  pdet ~ logleft2015 + participation + nep_mun, 
  data = cleaned_data, 
  method = "cem")
# summarize the matching results
# Appendix Tables 22 and 23
summary(cem_match2)
# extract balance data from the matchit object
balance_data2 <- summary(cem_match2)$sum.matched
# add subclass back to the cleaned dataset
cleaned_data$cem2_subclass <- cem_match2$subclass
# create a dataframe with only codigo and cem2_subclass from cleaned_data
cem_subclass_data <- cleaned_data[, c("codigo", "cem2_subclass")]
# merge cem2_subclass into data, keeping all rows from data
merged_data3 <- merge(
  data, 
  cem_subclass_data, 
  by = "codigo", 
  all.x = TRUE)
# observations in data but not in cleaned_data will automatically have NA in cem2_subclass.
str(merged_data3$cem2_subclass)
data <- merged_data3
```

```{r MergeData, message=TRUE, warning=FALSE}
# reshape from wide to long format
long_data <- data %>%
  pivot_longer(
    cols = c(
      starts_with("coca"), 
      starts_with("population"),
      starts_with("gdp"),
      starts_with("national_transfer"),
      starts_with("bombing"),
      starts_with("homicide")),
    names_to = c("variable", "year"),
    names_pattern = "(.*)(\\d{4})",
    values_to = "value"
  ) %>%
  pivot_wider(
    names_from = variable,
    values_from = value
  ) %>%
  mutate(year = as.integer(year))
# merge with annual assassinated data 
data6 <- merge(long_data, data5, by = c("codigo", "year"), all.x = TRUE)
# replace NA with 0: municipalities that were in data but were not in data0 did not have any assassinations
data7 <- replace(data6, is.na(data6), 0)
sum(data7$leader_killed) # should be 1379
length(unique(data7$codigo)) # should be 1121
data7 <- data7 %>% mutate(post = if_else(year >= 2017, 1, 0)) 
# adding year_factor and assassination as binary
data7 <- data7 %>% mutate(year_factor = as.factor(year)) %>% 
  mutate(year_factor = relevel(year_factor, ref = "2016")) %>%
  mutate(leader_killed_binary = if_else(leader_killed > 0, 1, 0)) 
nrow(data7) # should be 7847
# make numeric
data7 <- data7 %>%
  mutate(across(c(coca, homicide, gdp, national_transfer, bombing), as.numeric))
# filter out non-contested municipalities
data7contested <- data7[data7$contested > 0,]
nrow(data7contested) # should be 1344
# filtering out non-pdet provinces
data7pdet <- data7 %>% filter(department %in% (data7 %>% filter(pdet == 1) %>% distinct(department) %>% pull(department)))
nrow(data7pdet) # should be 4676
# only keep pdet & adjacent pdet municipalities
data7adjacent <- data7 %>% filter(pdet_adjacent <= 1)
nrow(data7adjacent) # should be 1960 = 7 years * 280 municipalities are either pdet or adjacent to it.
```

```{r DescriptiveStats}
# descriptive statistics 
# Appendix Table 4
table1(~ leader_killed + population + gdp + poverty + national_transfer + inst_development + inst_performance + bombing + homicide + coca + factor(conflict_intensity) + factor(contested) + project_num + ocadpaz + pdet_funds + total_investment +  factor(pnis) + families_pnis + logleft2015 + participation + nep_mun | factor(pdet), data=data7)
# correlations table
# Appendix Table 5
selected_data <- data7[, c("leader_killed", "population", "gdp", "poverty", "national_transfer", "inst_development", "inst_performance", "bombing", "homicide", "coca", "conflict_intensity", "contested")]
correlation_matrix(selected_data, digits = 2, use="lower", replace_diagonal = T)
save_correlation_matrix(df=selected_data, filename = "corr.csv", digits=2, use="lower")
```

```{r DiDfixest}
# DiD with the fixest package
DiD0 <- feglm(
  leader_killed ~ pdet * post | codigo + year, cluster = "codigo", family = poisson(), data = data7)
summary(DiD0) # very similar results to glmmTMB
# add province FE
DiD1 <- feglm(
  leader_killed ~ pdet * post | department + codigo + year, cluster = "codigo", family = poisson(), data = data7)
summary(DiD1) # very similar results to glmmTMB
# add covariates
DiD2 <- feglm(leader_killed ~ pdet * post + log(coca + 1) + log(gdp + 1) + log(national_transfer + 1) + bombing + log(population) | department + codigo + year, cluster="codigo", family = poisson(), data=data7)
summary(DiD2) # very similar results to glmmTMB
# fenegbin does not converge because of 0s
DiD2 <- fenegbin(leader_killed ~ pdet * post | codigo + year, cluster="codigo", data=data7)
summary(DiD2)
# I use the glmmTMB package in the next section. Results are very similar to fixest, but glmmTMB accounts for overdispersion. 
```

```{r negativebinomialmodels}
# check the distribution of the outcome variable
# Appendix Figure 1
Ord_plot(data2$leader_killed)
goodfittest <- goodfit(data2$leader_killed, type = "nbinomial")
summary(goodfittest)
goodfittest2 <- goodfit(data2$leader_killed, type ="poisson")
summary(goodfittest2)
# Appendix Figure 2
distplot(data2$leader_killed, type = "nbinomial")
# Appendix Figure 3
distplot(data2$leader_killed, type = "poisson")
# without province fixed effects 
fit_nb0 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year),
  family = nbinom2,
  data = data7)
# test dispersion
testDispersion(fit_nb0)
simulationOutput <- simulateResiduals(fit_nb0)
testZeroInflation(simulationOutput)
# make predictions
predictions <- predict(fit_nb0, type = "response")
# observed values
observed <- data7$leader_killed
# calculate RMSE
rmse <- sqrt(mean((predictions - observed)^2))
print(rmse)
# fit the Zero-Inflated Negative Binomial Model with fixed effects
fit_zinb0 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year),
  ziformula = ~ 1,
  family = nbinom2,
  data = data7)
summary(fit_zinb0)
# Fit poisson
fit_poisson0 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year),
  ziformula = ~ 1,
  family = poisson(link = "log"),
  data = data7)
summary(fit_poisson0)
# Appendix Table 7
lrtest(fit_nb0, fit_zinb0) 
# adding the term for zero-inflation (one extra parameter) does not significantly improve the fit of the model.
# overdispersion test for the Poisson model
# fit a Negative Binomial model with province fixed effects
fit_nb1 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# add control variables
fit_nb2 <- glmmTMB(leader_killed ~ pdet * post + (1 | codigo) + (1 | year) + (1 | department) 
                   + log(population) + log(gdp + 1) + log(national_transfer + 1) + log(coca + 1) + bombing,
  family = nbinom2,
  data = data7)
# check multicollinearity
check_collinearity(fit_nb2)
# check models with and without covariates
lrtest(fit_nb1, fit_nb2) # model with covariates fits the data better
# check against the zero-inflated model
fit_zinb <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year) + (1 | department),
  ziformula = ~ 1,
  family = nbinom2,
  data = data7)
summary(fit_zinb)
# zero-inflated model with control variables
fit_zinb3 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year) + (1 | department) 
  + log(coca + 1) + log(gdp + 1) + log(national_transfer + 1) 
  + bombing + log(population),
  ziformula = ~ 1,
  family = nbinom2,
  data = data7)
summary(fit_zinb3)
# check the residuals
residuals_nb1 <- residuals(fit_nb1, type = "pearson")
fitted_values_nb1 <- fitted(fit_nb1)
# plot residuals vs fitted values
plot(fitted_values_nb1, residuals_nb1,
     main = "Residuals vs Fitted Values",
     xlab = "Fitted Values",
     ylab = "Residuals")
  abline(h = 0, col = "red")
# histogram of residuals
hist(residuals_nb1, breaks = 50, main = "Histogram of Residuals", xlab = "Residuals")
# check the residuals for the zero-inflated model
residuals_zinb3 <- residuals(fit_zinb3, type = "pearson")
fitted_values_zinb3 <- fitted(fit_zinb3)
# plot residuals vs fitted values
plot(fitted_values_zinb3, residuals_zinb3,
     main = "Residuals vs Fitted Values",
     xlab = "Fitted Values",
     ylab = "Residuals")
abline(h = 0, col = "red")
# compare zinb and nb, both with control variables
lrtest(fit_nb1, fit_zinb3)
# check against the Poisson model
fit_poisson <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year) + (1 | department),
  family = poisson(link = "log"),
  data = data7)
summary(fit_poisson)
# Likelihood Ratio Test
lrtest(fit_nb1, fit_poisson) 
# The significant p-value (< 2.2e-16) indicates that the negative binomial model (fit_nb) provides a significantly better fit to the data compared to the Poisson model (fit_poisson), likely because it better accounts for overdispersion (the variance being greater than the mean) in the count outcome.
lrtest(fit_nb1, fit_zinb)
# The chi-squared value of 0 indicates that the addition of the zero-inflation component does not improve the model fit. The p-value of 0.9995, greater than the typical significance level of 0.05. This means there is no significant difference between the two models.
# Residual Analysis for the standard Negative Binomial model
residuals_nb <- residuals(fit_nb1)
# plot residuals
# Appendix Figure 4
plot(residuals_nb1, main = "Residuals of the NB Model", ylab = "Residuals", xlab = "Index")
abline(h = 0, col = "red")
# The red horizontal line at 0 represents where residuals should ideally be centered if the model fits the data well. Residuals should ideally be randomly scattered around this line with no clear pattern. 
# Plot residuals
plot(residuals_zinb3, main = "Residuals of the Zinb Model", ylab = "Residuals", xlab = "Index")
abline(h = 0, col = "red")
# check for overdispersion
# extract the residual deviance and degrees of freedom
residual_deviance <- sum(residuals(fit_nb1, type = "pearson")^2)
degrees_of_freedom <- df.residual(fit_nb1)
# calculate the dispersion parameter
dispersion_parameter <- residual_deviance / degrees_of_freedom
# print the results
cat("Residual Deviance:", residual_deviance, "\n")
cat("Degrees of Freedom:", degrees_of_freedom, "\n")
cat("Dispersion Parameter:", dispersion_parameter, "\n")
# check if the dispersion parameter indicates overdispersion
if (dispersion_parameter > 1) {
  cat("The model shows evidence of overdispersion.\n")
} else {
  cat("The model does not show evidence of overdispersion.\n")
} # the model does not show evidence of overdispersion. 
# only keeping provinces with pdet in data
fit_nb3 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7pdet)
# compare pdet and contiguous-to-pdet municipalities only
fit_nb4 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7adjacent)
# define the function to calculate RMSE
calculate_rmse <- function(model, data, response_var) {
  # make predictions
  predictions <- predict(model, type = "response")
  # observed values
  observed <- data[[response_var]]
  # calculate RMSE
  rmse <- sqrt(mean((predictions - observed)^2))
  return(rmse)
}
rmse_nb0 <- calculate_rmse(fit_nb0, data7, "leader_killed")
rmse_nb1 <- calculate_rmse(fit_nb1, data7, "leader_killed")
rmse_nb2 <- calculate_rmse(fit_nb2, data7, "leader_killed")
rmse_nb3 <- calculate_rmse(fit_nb3, data7, "leader_killed")
rmse_nb4 <- calculate_rmse(fit_nb4, data7, "leader_killed")
# print RMSE values
print(rmse_nb0)
print(rmse_nb1)
print(rmse_nb2)
print(rmse_nb3)
print(rmse_nb4)
# Manuscript Table 1 and Appendix Table 6
texreg(list(fit_nb0, fit_nb1, fit_nb2, fit_nb3, fit_nb4))
```

```{r PreTreatmentEffects}
# test parallel trends
fit_nb5 <- glmmTMB(
  leader_killed ~ pdet * year_factor + (1 | codigo),
  family = nbinom2,
  data = data7)
# test parallel trends
fit_nb6 <- glmmTMB(
  leader_killed ~ pdet * year_factor + (1 | codigo) + (1 | department),
  family = nbinom2,
  data = data7)
# add covariates
fit_nb7 <- glmmTMB(
  leader_killed ~ pdet * year_factor + log(coca + 1) + log(gdp + 1) + 
    log(national_transfer + 1) + bombing + log(population) + 
    (1 | codigo) + (1 | department),
  family = nbinom2,
  data = data7)
# pdet provinces only
fit_nb8 <- glmmTMB(
  leader_killed ~ pdet * year_factor + (1 | codigo) + (1 | department),
  family = nbinom2,
  data = data7pdet) 
# compare with adjacent municipalities
fit_nb9 <- glmmTMB(
  leader_killed ~ pdet * year_factor + (1 | codigo) + (1 | department),
  family = nbinom2,
  data = data7adjacent)
summary(fit_nb10)
# Appendix Table 8
texreg(list(fit_nb5, fit_nb6, fit_nb7, fit_nb8, fit_nb9))
```

```{r Logit}
# fit the logit model
fit_logit0 <- feglm(leader_killed_binary ~ pdet * post | codigo + year, cluster="codigo", 
      family = binomial(link = "logit"), 
      data = data7)
summary(fit_logit0)
# with province fixed effects
fit_logit1 <- feglm(leader_killed_binary ~ pdet * post | department + codigo + year, cluster="codigo", 
      family = binomial(link = "logit"), 
      data = data7)
summary(fit_logit1)
# with controls
fit_logit2 <- feglm(leader_killed_binary ~ pdet * post + log(coca + 1) + log(gdp + 1) + log(national_transfer + 1) + 
  bombing + log(population) | codigo + year + department, cluster="codigo", 
      family = binomial(link = "logit"), 
      data = data7)
summary(fit_logit2)
# Appendix Table 9
texreg(list(fit_logit0, fit_logit1, fit_logit2))
```

```{r Placebo}
# assign placebo treatment (before actual treatment period)
# only keep pre-treatment years
data7_filtered <- data7 %>%
  filter(year >= 2014 & year <= 2016) 
# assign 2014 as treatment year
data7_placebo <- data7_filtered %>%
  mutate(placebo = if_else(year > 2014, 1, 0))
# run model
fit_nb_placebo <- glmmTMB(
  leader_killed ~ pdet * placebo + (1 | codigo) + (1 | department),
  family = nbinom2,
  data = data7_placebo
)
summary(fit_nb_placebo)
# assign 2015 as treatment year
data7_placebo2 <- data7_filtered %>%
  mutate(placebo = if_else(year > 2015, 1, 0))
# run model
fit_nb_placebo2 <- glmmTMB(
  leader_killed ~ pdet * placebo + (1 | codigo) + (1 | department),
  family = nbinom2,
  data = data7_placebo2
)
summary(fit_nb_placebo2)
# Appendix Table 10
texreg(list(fit_nb_placebo, fit_nb_placebo2))
```

```{r SensitivityAnalysis}
# sensitivity analysis
# remotes::install_github("asheshrambachan/HonestDiD")
# event study 
event_study <- feols(leader_killed ~ i(year, pdet, ref= "2016") | codigo + year , cluster="codigo", data=data7adjacent)
# Appendix Figure 5
iplot(event_study)
etable(event_study)
betahat <- summary(event_study)$coefficients
sigma <- summary(event_study)$cov.scaled
# event study with covariates
event_study2 <- feols(leader_killed ~ i(year, pdet, ref= "2016") + log(coca + 1) + log(gdp + 1) + log(national_transfer + 1) + bombing + log(population)  | codigo + year + department^year, cluster="codigo", data=data7)
# Manuscript Figure 3
iplot(event_study2)
# extract the summary of the model
summary_event_study2 <- summary(event_study2)
# extract all coefficients
all_coefficients <- summary_event_study2$coefficients
interaction_terms <- c("year::2014:pdet", "year::2015:pdet", "year::2017:pdet", 
                       "year::2018:pdet", "year::2019:pdet", "year::2020:pdet")
betahat2 <- all_coefficients[interaction_terms]
sigma2 <- summary_event_study2$cov.scaled[interaction_terms, interaction_terms]
# sensitivity analysis for the average over post periods (adjacent municipalities)
didsens1 <- createSensitivityResults_relativeMagnitudes(
  betahat = betahat,
  sigma = sigma,
  numPrePeriods = 2,
  numPostPeriods = 4,
  Mbarvec = seq(0.5,2,by=0.5),
  l_vec = rep(0.2,4),
  grid.lb = -1,grid.ub=1,gridPoints=2000)
# get results
didsens1
# plot
es_results1 <- constructOriginalCS(betahat = betahat,
                                   sigma = sigma,
                                   numPrePeriods = 2,
                                   numPostPeriods = 4,
                                   l_vec = rep(0.2,4))

createSensitivityPlot_relativeMagnitudes(didsens1,es_results1)
# Sensitivity analysis for the average over post periods (all municipalities with covariates)
didsens2 <- createSensitivityResults_relativeMagnitudes(
  betahat = betahat2,
  sigma = sigma2,
  numPrePeriods = 2,
  numPostPeriods = 4,
  Mbarvec = seq(0.5,2,by=0.5),
  l_vec = rep(0.2,4),
  grid.lb = -1,grid.ub=1,gridPoints=2000)
# get results
didsens2
# plot
es_results2 <- constructOriginalCS(betahat = betahat2,
                                   sigma = sigma2,
                                   numPrePeriods = 2,
                                   numPostPeriods = 4,
                                   l_vec = rep(0.2,4))
# Manuscript Figure 4
createSensitivityPlot_relativeMagnitudes(didsens2,es_results2)
```

```{r HeterogenousEffects}
# Heterogeneous Treatment Effects
# Manuscript Table 4 and Appendix Table 12
# coca
fit_het1 <- glmmTMB(
  leader_killed ~  pdet * post * log(coca + 1) + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# poverty
fit_het2 <- glmmTMB(
  leader_killed ~  pdet * post * poverty + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# homicides
fit_het3 <- glmmTMB(
  leader_killed ~  pdet * post * log(homicide+1) + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# contested (binary) only significant at 0.1
fit_het4 <- glmmTMB(
  leader_killed ~  pdet * post * contested + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# institutional capacity: Gasto, national transfer, SGR not sig
fit_het5 <- glmmTMB(
  leader_killed ~  pdet * post * log(gdp + 1) + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# Manuscript Table 4 and Appendix Table 12
texreg(list(fit_het1, fit_het2, fit_het3, fit_het4, fit_het5))
```

```{r ContestedSample}
# In contested municipalities
fit_nbcon0 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year),
  family = nbinom2,
  data = data7contested)
# including province FE
fit_nbcon1 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7contested)
# including covariates
fit_nbcon2 <- glmmTMB(
  leader_killed ~ pdet * post + log(coca + 1) + log(gdp + 1) + log(national_transfer + 1) + 
  bombing + log(population) + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7contested)
# Manuscript Table 3 and Appendix Table 13
texreg(list(fit_nbcon0, fit_nbcon1, fit_nbcon2))
```

```{r Investment}
# effect of various programs and amount/sources of funding
# no. of projects
fit_invest0 <- glmmTMB(
leader_killed ~  post * project_num + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
summary(fit_invest0)
# total amount of investment
fit_invest1 <- glmmTMB(
leader_killed ~  post * log(total_investment + 1) + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# total pdet funds
fit_invest2 <- glmmTMB(
leader_killed ~  post * log(pdet_funds + 1) + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# ocadpaz funds
fit_invest3 <- glmmTMB(
leader_killed ~  post * log(ocadpaz+1) + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# families participating in the PNIS program
fit_invest4 <- glmmTMB(
leader_killed ~  post * log(families_pnis + 1) + (1 | codigo) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# Manuscript Table 2 and Appendix Table 14
texreg(list(fit_invest0, fit_invest1, fit_invest2, fit_invest3, fit_invest4))
```

```{r MatchedSamples}
# matched with propensity scores
fit_matched0 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | fm0) + (1 | year),
  family = nbinom2,
  data = data7)
# add province FE
fit_matched1 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | fm0) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# add covariates
fit_matched2 <- glmmTMB(
  leader_killed ~ pdet * post + log(gdp + 1) + log(national_transfer + 1) + 
  bombing + log(population) + (1 | fm0) + (1 | year) + (1 | department),
  family = nbinom2,
  data = data7)
# Appendix Table 18
texreg(list(fit_matched0, fit_matched1, fit_matched2))

# matched with cem: criminal pathway
fit_matched3 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | cem_subclass) + (1 | year) ,
  family = nbinom2,
  data = data7)
# Appendix Table 22
texreg(fit_matched3)
# matched with cem with Albarracin et al. (2023) variables: political pathway
fit_matched4 <- glmmTMB(
  leader_killed ~ pdet * post + (1 | cem2_subclass) + (1 | year),
  family = nbinom2,
  data = data7)
# Appendix Table 23
texreg(fit_matched4)
```

```{r Figure2}
# create a monthly dataset to get the year-half graph
# Create a data frame with all combinations of codigo, year, and month
codigo <- unique(data$codigo)
years <- 2014:2020
months <- 1:12
combinations <- expand.grid(codigo = codigo, year = years, month = months)
# Merge the combinations with the subsetted data
final_data <- merge(combinations, data, by = "codigo", all.x = TRUE)
# Merge with assassinated data (data2)
merged_data <- merge(final_data, data3, by = c("codigo", "year", "month"), all.x = TRUE)
length(unique(merged_data$codigo)) # should be 1121
sum(merged_data$leader_killed, na.rm = TRUE) # should be 1379
merged_data <- replace(merged_data, is.na(merged_data), 0)
# create year-half
merged_data$year_half <- ifelse(merged_data$month <= 6, paste0(merged_data$year, "-1"), paste0(merged_data$year, "-2"))
merged_data <- merged_data %>%
  group_by(year_half, pdet) %>% summarize(total_killed = sum(leader_killed, na.rm = TRUE)) %>% ungroup() 
# Manuscript FIgure 2
OO <- ggplot(merged_data, aes(x = year_half, y = total_killed, group = pdet, color = as.factor(pdet))) +
  geom_line() +
  geom_vline(xintercept = "2016-2", linetype = "dashed") +
  scale_color_manual(values = c("grey", "black")) +
  labs(
    title = "Selective Violence by PDET Status",
    x = "Half-year",
    y = "Assassinated leaders",
    color = "PDET"
  ) +
  theme_classic(base_family = "Times New Roman") + 
  theme(
    axis.line.x = element_line(color = "black", size = 0.5),
    axis.line.y = element_line(color = "black", size = 0.5),
    axis.text.x = element_text(color = "black", size = 12, face = "plain", family = "Times New Roman", angle = 90, hjust = 1),
    axis.text.y = element_text(color = "black", size = 12, face = "plain", family = "Times New Roman"),
    plot.title = element_text(family = "Times New Roman", face = "bold", size = 14),
    legend.title = element_text(family = "Times New Roman", face = "bold", size = 12),
    legend.text = element_text(family = "Times New Roman", size = 10)
  )
OO
file_path <- "OO.jpg" 
# Save as JPEG with good quality
ggsave(file_path, plot = OO, width = 8, height = 6, dpi = 1000)
```

```{r Map}
# Maps 
# devtools::install_github("nebulae-co/colmaps")
colnames(municipios@data)
# reading the shapefiles
sf <- st_read(dsn="COL_adm", layer="COL_adm2")
shape <- readOGR(dsn="COL_adm", layer="COL_adm2")
# plotting the shapefile
plot(shape)
plot(sf)
# plot the districts only
plot(sf["NAME_1"], axes = TRUE, main = "Districts")
unique(shape@data$NAME_1)
# Create a new column in the municipios data to identify pdet municipalities
municipios@data <- municipios@data %>%
  mutate(pdet = ifelse(id %in% pdet, "pdet", "Non-pdet"))
# Convert municipios to an sf object
municipios_sf <- st_as_sf(municipios)
# plot the map
ggplot(data = municipios_sf) +
  geom_sf(aes(fill = pdet), color = "black") +
  scale_fill_manual(values = c("pdet" = "grey", "Non-pdet" = "white")) +
  theme_void() +
  labs(title = "pdet Municipalities in Colombia")
# create a new column in the municipios data 
municipios@data$pdet <- ifelse(municipios@data$id %in% pdet, "pdet", "Non-pdet")
# determine provinces with at least one pdet municipality
provinces_with_pdet <- municipios@data %>%
  filter(pdet == "pdet") %>%
  pull(id_depto) %>%
  unique()
# ddd a column to identify if a province has at least one pdet municipality
municipios@data$province_has_pdet <- ifelse(municipios@data$id_depto %in% provinces_with_pdet, "yes", "no")
# convert municipios to an sf object 
municipios_sf <- st_as_sf(municipios)
# generate province borders
provinces_sf <- municipios_sf %>%
  group_by(id_depto) %>%
  summarise(geometry = st_union(geometry), .groups = 'drop')
# Manuscript Figure 5
# plot the map with conditional borders
allmap <- ggplot() +
  geom_sf(data = municipios_sf %>% filter(province_has_pdet == "yes"), aes(fill = pdet), color = "black", linewidth = 0.2) +  # Thin borders for municipalities in provinces with pdet
  geom_sf(data = municipios_sf %>% filter(province_has_pdet == "no"), aes(fill = pdet), color = NA) +  # No municipality borders for other provinces
  geom_sf(data = provinces_sf, fill = NA, color = "black", linewidth = 0.4) +  # Province borders
  scale_fill_manual(values = c("pdet" = "grey", "Non-pdet" = "white")) +
  theme_void() +
  labs(title = "pdet Municipalities in Colombia")
ggsave("pdet_Municipalities_in_Colombia.jpg", plot = allmap, width = 10, height = 10, dpi = 300)
```

```{r MapNarino}
# Narino
nariño_sf <- sf %>%
  filter(NAME_1 == "Nariño")
# in color
plot(nariño_sf["NAME_2"], axes = TRUE, main = "Nariño - Districts")
# in grey
plot(nariño_sf["NAME_2"], col = "lightgrey", axes = TRUE, main = "Nariño - Municipalities")
nariño_sf$NAME_2
# add pdet
# list of municipalities to be marked as 1
pdet <- c("Cumbitara", "El Rosario", "Leiva", "Los Andes", "Policarpa",
                           "Barbacoas", "El Charco", "La Tola", "Magüí", "Mosquera", 
                           "Olaya Herrera", "Francisco Pizarro", "Ricaurte", "Roberto Payán", "Santa Bárbara", "Tumaco")
# create the new column
nariño_sf$pdet <- ifelse(nariño_sf$NAME_2 %in% pdet, 1, 0)
# filter the data for pdet == 1
pdet_sf <- nariño_sf[nariño_sf$pdet == 1, ]
# calculate a point on the surface for the filtered municipalities
pdet_points <- st_point_on_surface(pdet_sf)
# create a data frame for labeling with numbers
label_data <- data.frame(
  x = st_coordinates(pdet_points)[, 1],
  y = st_coordinates(pdet_points)[, 2],
  label = 1:nrow(pdet_sf))
# plot the map with numbers on municipalities
QQ <- ggplot() +
  geom_sf(data = nariño_sf, aes(fill = ifelse(pdet == 1, "pdet", "Not pdet")), color = "black", size = 0.5) +
  geom_text(data = label_data, aes(x = x, y = y, label = label), color = "black", size = 3) +
  scale_fill_manual(values = c("pdet" = "lightgrey", "Not pdet" = "white"),
                    labels = c("pdet" = "pdet", "Not pdet" = "Non-pdet"),
                    guide = guide_legend(title = "")) +
  labs(title = "Nariño") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78))
file_path <- "QQ.jpg" 
# Save as JPEG with good quality
ggsave(file_path, plot = QQ, width = 6, height = 6, dpi = 1000)
## graph municipalities with assassination counts
# change the municipality names to match assassination data
changes <- tribble(
  ~old_name, ~new_name,
  "Ancuyá", "Ancuya",
  "Cuaspud", "Cuaspud Carlosama",
  "El Tablón de Gomez", "El Tablón De Gómez",
  "La Unión de Sucre", "La Unión",
  "San Juan de Pasto", "Pasto",
  "San Pedro de Cartago", "San Pedro De Cartago",
  "Santa Cruz", "Santacruz",
  "Tumaco", "San Andrés De Tumaco"
)
# apply the changes to nariño_sf
nariño_sf <- nariño_sf %>%
  left_join(changes, by = c("NAME_2" = "old_name")) %>%
  mutate(NAME_2 = if_else(is.na(new_name), NAME_2, new_name)) %>%
  select(-new_name)
# merge with assassinations data
NarinoKilled <- data6 %>%
  filter(department == "Narino")
# count the total number of leader_killed for each municipality
# wartime numbers
NarinoKilled3 <- NarinoKilled %>%
  group_by(codigo) %>% filter(year < 2017) %>%
  summarise(total_killed = sum(leader_killed, na.rm = TRUE)) %>% ungroup()
# merge
nariño_sf1 <- nariño_sf %>%
  left_join(NarinoKilled3, by = c("NAME_2" = "municipality"))
### create a map with leader killed
# generate a gradient color palette with 99 shades of grey
grey_palette <- c("white", rev(scales::grey_pal()(50)))
# define two separate grey palettes
grey_palette_1 <- c("white", rev(grey_pal()(25)[1:12])) # The darker half of the grey scale
grey_palette_2 <- c("white", rev(grey_pal()(25)[13:25])) # The lighter half of the grey scale
# create the plot 
PQ <- ggplot() +
  geom_sf(data = nariño_sf1, aes(fill = total_killed)) +
  scale_fill_gradientn(colors = grey_palette_2, 
  name = "Total Killed",  guide = guide_colorbar(title.position = "top", title.hjust = 0.5, label = "")) +
  labs(title = "Assassinations during Conflict") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78))   
PQ
file_path1 <- "PQ.jpg" 
ggsave(file_path1, plot = PQ, width = 6, height = 6, dpi = 1000)
# postwar graph
NarinoKilled4 <- NarinoKilled %>%
  group_by(codigo) %>% filter(year > 2016) %>%
  summarise(total_killed = sum(leader_killed, na.rm = TRUE))  %>% ungroup()
# merge
nariño_sf2 <- nariño_sf %>%
  left_join(NarinoKilled4, by = c("NAME_2" = "municipality"))
### create a map with leader killed
PP <- ggplot() +
  geom_sf(data = nariño_sf2, aes(fill = total_killed)) + 
  geom_sf(data = nariño_sf2 %>% filter(pdet == 1), fill = NA, color = "black", linewidth = 0.7) +
  scale_fill_gradientn(colors = grey_palette_1, 
  name = "Total Killed",  guide = guide_colorbar(title.position = "top", title.hjust = 0.5, label = "")) +
  labs(title = "Post-Conflict Assassinations") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78))  
PP
file_path <- "PP.jpg" 
ggsave(file_path, plot = PP, width = 6, height = 6, dpi = 1000)
# combine the wartime and post-war plots
# Manuscript Figure 6 - Nariño
grid.arrange(PQ, PP, ncol = 2)
combined_plot <- grid.arrange(PQ, PP, ncol = 2)
file_path <- "combined_plot.jpg" 
ggsave(file_path, plot = combined_plot, width = 10, height = 6, dpi = 1000)
```

```{r MapAntioquia}
# Antioquia
Antioquia_sf <- sf %>%
  filter(NAME_1 == "Antioquia")
Antioquia_sf$NAME_2
# add pdet
# list of municipalities to be marked with 1
pdet_adjacent <- c("Amalfi", "Anorí", "Briceño", "Cáceres", "Caucasia", "El Bagre", "Ituango", "Nechí", "Remedios", "Segovia", "Tarazá", "Valdivia", "Zaragoza", "Murindó", "Vigía del Fuerte", "Yondó", "Apartadó", "Carepa", "Chigorodó", "Dabeiba", "Mutatá", "Necoclí", "San Pedro de Urabá", "Turbo")
# create the new column
Antioquia_sf$pdet <- ifelse(Antioquia_sf$NAME_2 %in% pdet_adjacent, 1, 0)
sum(Antioquia_sf$pdet) # should be 24
# create the previous map with ggplot
# filter the data for pdet == 1
pdet_sf2 <- Antioquia_sf[Antioquia_sf$pdet == 1, ]
# calculate a point on the surface for the filtered municipalities
pdet_points2 <- st_point_on_surface(pdet_sf2)
# create a data frame for labeling with numbers
label_data2 <- data.frame(
  x = st_coordinates(pdet_points2)[, 1],
  y = st_coordinates(pdet_points2)[, 2],
  label = 1:nrow(pdet_sf2))
# Plot the map with numbers on municipalities
QQQ <- ggplot() +
  geom_sf(data = Antioquia_sf, aes(fill = ifelse(pdet == 1, "pdet", "Not pdet")), color = "black", size = 0.5) +
  geom_text(data = label_data2, aes(x = x, y = y, label = label), color = "black", size = 3) +
  scale_fill_manual(values = c("pdet" = "lightgrey", "Not pdet" = "white"),
                    labels = c("pdet" = "pdet", "Not pdet" = "Non-pdet"),
                    guide = guide_legend(title = "")) +
  labs(title = "Antioquia") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78))
file_path <- "QQQ.jpg" 
# save as JPEG with good quality
ggsave(file_path, plot = QQQ, width = 6, height = 6, dpi = 1000)
## graph municipalities with assassination counts
# merge with assassinations data
AntioquiaKilled <- data6 %>%
  filter(department == "Antioquia")
# change the municipality names to match assassination data
changes2 <- tribble(
  ~old_name, ~new_name,
  "Santa Fe de Antioquia", "Santa Fé De Antioquia",
  "San Andrés de Cuerquia", "San Andrés De Cuerquía",
  "Pequé", "Peque",
  "Carolina del Principe", "Carolina",
  "Don Matías", "Donmatías",
  "El Carmen de Viboral", "El Carmen De Viboral",
  "El Santuario", "El Santuario",
  "Vigía del Fuerte", "Vigía Del Fuerte",
  "La Unión de Sucre", "La Unión",
  "San Pedro de los Milagros", "San Pedro De Los Milagros",
  "San Pedro de Urabá", "San Pedro De Urabá",
  "San José de la Montaña", "San José De La Montaña",
  "San Juan de Urabá", "San Juan De Urabá")
# apply the changes to nariño_sf
Antioquia_sf <- Antioquia_sf %>%
  left_join(changes2, by = c("NAME_2" = "old_name")) %>%
  mutate(NAME_2 = if_else(is.na(new_name), NAME_2, new_name)) %>%
  select(-new_name)
# count the total number of leader_killed for each municipality
# wartime numbers
AntioquiaKilled3 <- AntioquiaKilled %>%
  group_by(codigo) %>% filter(year < 2017) %>%
  summarise(total_killed = sum(leader_killed, na.rm = TRUE))  %>% ungroup()
# merge
Antioquia_sf1 <- Antioquia_sf %>%
  left_join(AntioquiaKilled3, by = c("NAME_2" = "municipality"))
### create a map with leader killed
# generate a gradient color palette with 99 shades of grey
grey_palette <- c("white", rev(scales::grey_pal()(50)))
library(scales)
# define two separate grey palettes
grey_palette_1 <- c("white", rev(grey_pal()(25)[1:12])) # The darker half of the grey scale
grey_palette_2 <- c("white", rev(grey_pal()(25)[13:25])) # The lighter half of the grey scale
# create the plot 
PQQ <- ggplot() +
  geom_sf(data = Antioquia_sf1, aes(fill = total_killed)) +
  scale_fill_gradientn(colors = grey_palette_2, 
  name = "Total Killed",  guide = guide_colorbar(title.position = "top", title.hjust = 0.5, label = "")) +
  labs(title = "Assassinations during Conflict") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78))  
PQQ
file_path1 <- "PQQ.jpg" 
ggsave(file_path1, plot = PQ, width = 6, height = 6, dpi = 1000)
# postwar graph
AntioquiaKilled4 <- AntioquiaKilled %>%
  group_by(codigo) %>% filter(year > 2016) %>%
  summarise(total_killed = sum(leader_killed, na.rm = TRUE)) %>% ungroup()
# merge
Antioquia_sf2 <- Antioquia_sf %>%
  left_join(AntioquiaKilled4, by = c("NAME_2" = "municipality"))
### create a map with leader killed
PPQ <- ggplot() +
  geom_sf(data = Antioquia_sf2, aes(fill = total_killed)) +
  geom_sf(data = Antioquia_sf2 %>% filter(pdet == 1), fill = NA, color = "black", linewidth = 0.7) +
  scale_fill_gradientn(colors = grey_palette_1, 
  name = "Total Killed", guide = guide_colorbar(title.position = "top", title.hjust = 0.5, label = "")) +
  labs(title = "Post-Conflict Assassinations") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78)) 
PPQ 
file_path <- "/PPQ.jpg" 
ggsave(file_path, plot = PP, width = 6, height = 6, dpi = 1000)
# combine the wartime and post-war plots
# Manuscript Figure 6 - Antioquia 
grid.arrange(PQQ, PPQ, ncol = 2)
combined_plot2 <- grid.arrange(PQQ,PPQ, ncol = 2)
file_path <- "combined_plot2.jpg" 
ggsave(file_path, plot = combined_plot2, width = 10, height = 6, dpi = 1000)
```

```{r MapCatatumbo}
# Norte de Santander
Catatumbo_sf <- sf %>%
  filter(NAME_1 == "Norte de Santander")
Catatumbo_sf$NAME_2
# add pdet
# list of municipalities to be marked with 1
pdet3 <- c("Convención", "El Carmen", "El Tarra", "Hacarí", "San Calixto", "Sardinata", "Teorama", "Tibú")
# creating the new column
Catatumbo_sf$pdet <- ifelse(Catatumbo_sf$NAME_2 %in% pdet3, 1, 0)
sum(Catatumbo_sf$pdet) # should be 24
# filter the data for pdet == 1
pdet_sf3 <- Catatumbo_sf[Catatumbo_sf$pdet == 1, ]
# calculate a point on the surface for the filtered municipalities
pdet_points3 <- st_point_on_surface(pdet_sf3)
# create a data frame for labeling with numbers
label_data3 <- data.frame(
  x = st_coordinates(pdet_points3)[, 1],
  y = st_coordinates(pdet_points3)[, 2],
  label = 1:nrow(pdet_sf3))
# plot the map with numbers on municipalities
QQQQ <- ggplot() +
  geom_sf(data = Catatumbo_sf, aes(fill = ifelse(pdet == 1, "pdet", "Not pdet")), color = "black", size = 0.5) +
  geom_text(data = label_data3, aes(x = x, y = y, label = label), color = "black", size = 3) +
  scale_fill_manual(values = c("pdet" = "lightgrey", "Not pdet" = "white"),
                    labels = c("pdet" = "pdet", "Not pdet" = "Non-pdet"),
                    guide = guide_legend(title = "")) +
  labs(title = "Norte de Santander") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78))

file_path <- "QQQQ.jpg" 
# Save as JPEG with good quality
ggsave(file_path, plot = QQQQ, width = 6, height = 6, dpi = 1000)
## Graph municipalities with assassination counts
# Merge with assassinations data
CatatumboKilled <- data6 %>%
  filter(department == "Norte Santander")
# change the municipality names to match assassination data
changes3 <- tribble(
  ~old_name, ~new_name,
  "Abrego", "Ábrego",
  "La Playa de Belén", "La Playa",
  "Salazar de las Palmas", "Salazar",
  "San José de Cúcuta", "San José De Cúcuta",
  "Santo Domingo de Silos", "Silos",
  "Villa del Rosario", "Villa Del Rosario")
# apply the changes to nariño_sf
Catatumbo_sf <- Catatumbo_sf %>%
  left_join(changes3, by = c("NAME_2" = "old_name")) %>%
  mutate(NAME_2 = if_else(is.na(new_name), NAME_2, new_name)) %>%
  select(-new_name)
# count the total number of leader_killed for each municipality
# wartime numbers
CatatumboKilled3 <- CatatumboKilled %>%
  group_by(codigo) %>% filter(year < 2017) %>%
  summarise(total_killed = sum(leader_killed, na.rm = TRUE))  %>% ungroup()
# merge
Catatumbo_sf1 <- Catatumbo_sf %>%
  left_join(CatatumboKilled3, by = c("NAME_2" = "municipality"))
### create a map with leader killed
grey_palette <- c("white", rev(scales::grey_pal()(50)))
grey_palette_1 <- c("white", rev(grey_pal()(25)[1:12])) 
grey_palette_2 <- c("white", rev(grey_pal()(25)[13:25])) 
# Create the plot 
PQQQ <- ggplot() +
  geom_sf(data = Catatumbo_sf1, aes(fill = total_killed)) +
  scale_fill_gradientn(colors = grey_palette_2, 
  name = "Total Killed",  guide = guide_colorbar(title.position = "top", title.hjust = 0.5, label = "")) +
  labs(title = "Assassinations during Conflict") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78))  
PQQQ
file_path1 <- "PQQQ.jpg" 
ggsave(file_path1, plot = PQQQ, width = 6, height = 6, dpi = 1000)
# postwar graph
CatatumboKilled4 <- CatatumboKilled %>%
  group_by(codigo) %>% filter(year > 2016) %>%
  summarise(total_killed = sum(leader_killed, na.rm = TRUE))  %>% ungroup()
# merge
Catatumbo_sf2 <- Catatumbo_sf %>%
  left_join(CatatumboKilled4, by = c("NAME_2" = "municipality"))
### create a map with leader killed
PPQQ <- ggplot() +
  geom_sf(data = Catatumbo_sf2, aes(fill = total_killed)) +
  geom_sf(data = Catatumbo_sf2 %>% filter(pdet == 1), fill = NA, color = "black", linewidth = 0.7) +
  scale_fill_gradientn(colors = grey_palette_1, 
  name = "Total Killed", guide = guide_colorbar(title.position = "top", title.hjust = 0.5, label = "")) +
  labs(title = "Post-Conflict Assassinations") +
  theme_void() +
  theme(legend.justification = "left",  
        legend.position = c(0.78, 0.78)) 
PPQQ 
file_path <- "PPQ.jpg" 
ggsave(file_path, plot = PP, width = 6, height = 6, dpi = 1000)
# combine the wartime and post-war plots
# Manuscript Figure 6 - Norte de Santander
grid.arrange(PQQQ, PPQQ, ncol = 2)
combined_plot3 <- grid.arrange(PQQQ,PPQQ, ncol = 2)
file_path <- "combined_plot3.jpg" 
ggsave(file_path, plot = combined_plot3, width = 10, height = 6, dpi = 1000)
```